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ABSTRACT 

We use new deep near-infrared (NIR) and mid-infrared (MIR) observations to analyze the 850 /rm 
image of the Great Observatories Origins Deep Survey-North region around the Hubble Deep Field- 
North. We show that much of the submillimeter background at this wavelength is picked out by sources 
with H{AB) or 3.6 /rm(AH) < 23.25 (1.8 /iJy). These sources contribute an 850 ^m background of 
24 ± 2 Jy deg“^. This is a much higher fraction of the measured background (31 — 45 Jy deg“^) 
than is found with current 20 cm or 24 ^m samples. Roughly one-half of these NIR-selected sources 
have spectroscopic identifications, and we can assign robust photometric redshifts to nearly all of 
the remaining sources using their UV to MIR spectral energy distributions. We use the redshift 
and spectral type information to show that a large fraction of the 850 /rm background light comes 
from sources with z = 0 — 1.5 and that the sources responsible have intermediate spectral types. 
Neither the elliptical galaxies, which have no star formation, nor the bluest galaxies, which have little 
dust, contribute a significant amount of 850 /rm light, despite the fact that together they comprise 
approximately half of the galaxies in the sample. The galaxies with intermediate spectral types have 
a mean flux of 0.40 ± 0.03 mJy at 850 /im and 9.1 ± 0.3 /rJy at 20 cm. 

The redshift distribution of the NIR-selected 850 /rm light lies well below that of the much smaller 
amount of light traced by the more luminous, radio-selected submillimeter sources. We therefore 
require a revised star-formation history with a lower star-formation rate at high redshifts. We use a 
stacking analysis of the 20 cm light in the NIR sample to show that the star-formation history of the 
total 850 p,m sample is relatively flat down to z ~ 1 and that half of the total star formation occurs 
at redshifts z < 1.4. 

Subject headings: cosmology: observations — galaxies: evolution — galaxies: formation — galaxies : 
starburst — infrared: galaxies — submillimeter 


1. INTRODUCTION 

The integrated extragalactic background light (EBL) 
is a measure of the history of the luminous energy pro¬ 
duction of the universe from both star formation and 
active galactic nuclei (AGNs). Directly emitted light is 
seen in the UV and optical, whereas dust reradiated en¬ 
ergy appears in the far-infrared (FIR) and submillime¬ 
ter. COBE obtained detailed measurement s of the EBL 
at FIR and submillimet er wavelengths (e.g.. lPuget et al.l 
ll996tlFixsen et al.ll998^ . showing that the total radiated 
emission reprocessed by dust in the FIR/submillimeter is 
comparable to the total measured optical EBL. However, 
to proceed further, we also need to know the redshift dis¬ 
tribution of the sources contributing to the submillime¬ 
ter background, and this information has been extremely 
difficult to obtain. 

In the last decade, the submillimeter/millimeter EBL 
has been resolved into discrete sources by deep sur¬ 
veys with the Submillimeter Common-User Bolometer 
Array (SCUBA) on the 15 m James Clerk Maxwell 
Telescope (JCMT) and with the Max-Plank Millime¬ 
ter Bolometer array on the 30 m IRAM telescope. 
Blank-field surveys have resolved sources in the 2 — 
20 mJy range that account for ~ 20 — 30% of the 
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With the help of strong lensing, surveys in cluster 
fields have resolved sources over the 0.3 — 2 mJy range 
that account for a further 45 — 65% of the 850 um 
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EBL l|Smail, Ivison, fc Blain|199'^ [Chapman et ak l2002t 
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gether these surveys provide a nearly complete resolu¬ 
tion of the background at 850 ^m. The “typical” source 
contributing to the 850 fim EBL has a mean flux of 
about 0.9 mJv and a median flux of about 0.6 mJy 
llCowie. Barger, fc Kne.ibll20?il . 

However, the redshift follow-up of the submillime¬ 
ter sources has been very slow. Because of the large 
beam size (15") of SCUBA and the optically-faint 
nature of the dusty sources, identifying the optical 
and near-infrared (NIR) counterparts to the submil¬ 
limete r sources is time c onsuming (e.g., IBarger et a.U 
119994 E vison et al ] l2flflFll . To date, the most suc¬ 
cessful identifications of the submillimeter sources rely 
on the empirical correlation between the nonthermal 
radio emission and the thermal dust emission (e.g., 
iCondoiJ 1199211 . Once the radio counterparts to the 
submillimeter sources are detected by radio interfer¬ 
ometers, the redshifts of the sources can be crudely 
estimated using the radio-to-submillimeter flux ratios 
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| 2003bD or accura t ely me a sured with optical spectroscopy 
l|Cha,nma,n et ID [2M HID). The radio-identified 
sources are mostly bright (3> 2 mJy) submillimeter 
sources at z = 1.5 — 3.5, with properties similar 
to the local ultraluminous infrared galaxies (ULIRGs; 
Lir, > 10^^ Lq, where Lt r is the 8 — 1000 »m in - 
frared luminosity; see, e.g., ISanders fc Mirab’^IlQQfiH . 
We note, however, that because of the RT-correction 
and the sensitivity limit in the radio, only 60% of 
the bright submillimeter sources are i dentified in the ra¬ 
dio iBarger. Cowie. fc RichardsI 1200(11 . It is not known 
whether the remaining 40% are at higher redshifts that 
simply cannot be reached by current radio telescopes. 

Importantly, however, the properties and redshift dis¬ 
tribution of the faint submillimeter sources that domi¬ 
nate the submillimeter EBL remain essentially unknown. 
The absence of any redshift information for more than 
90% of the 850 p,m EBL represents a formidable uncer¬ 
tainty in determining the star-formation history, and this 
is what we aim to resolve in the present paper. 

Like the radio emission, the mid-infrared (MIR) emis¬ 
sion at > 5 /im could serve as another proxy to 

the submillimeter emission, since it also comes from 
dust. The MIR window has been opened by th e In¬ 
frared Array Camera (IRAC, iFazio et~a D l2(^ and 
the Multiband Im aging Photometer for Spitzer (MIPS, 
Rieke et al.ll2004ll on the Spitzer Space T el escoye (e.g., 
Huang et al.1 12004 iSerieant et al.f 12004 llvison et alJ 
20f)4nEgami~er^,l J i2004|i ! MIPS should be sensitive to 
z < I galaxies with infrared luminosities similar to lo¬ 
cal normal galaxies (Lir ~ 10^° Lq, corresponding to 
~ 0.1 mJy at 850 /rm) and to z < 3.5 ULIRGs (i.e., typ¬ 
ical of the bright submillimeter sources). Thus, MIPS 
should be able to detect the radio-identified submillime¬ 
ter sources at z < 3.5 and to provide a large sample of 
faint sources that are beyond the co nfusion limit of cur - 
rent submillimeter telescopes (e.g.. IcTha.rv et ;^l2flfl4jl . 
However, as we shall show in this paper, even the extraor¬ 
dinarily deep MIPS data of the Great Observatories Ori¬ 
gins Deep Survey-North (GOODS-N) Spitzer Legacy Sci¬ 
ence Program in the Hubble Deep Field-North (HDF-N) 
region does not substantially identify the 850 /rm EBL. 

Remarkably, however, the combination of a J or H- 
band sample (selected from images obtained with the 
new generation of ground-based, wide-field NIR cam¬ 
eras) and the IRAC 3.6 /im sample does identify much 
of the 850 /im EBL. We show this using t he iL-band im¬ 
age of the GOODS-N region obtained bv iTrouille et al.1 
ll2006tl . This result makes sense if the bulk of the sources 
contributing to the 850 /im EBL are actually at lower 
redshifts and luminosities than those identified at the 
brighter submillimeter fluxes. Such sources have strong 
rest-frame optical/NIR counterparts that are picked up 
in the NIR sample. We utilize the spectroscopic and 
photometric redshift information on our NIR sample to 
confirm this result. We find that more than half of the 
850 /im EBL arises in sources with z < 1.5 and that the 
sources that are responsible have intermediate spectral 
types. Neither the elliptical galaxies, which have no star 
formation, nor the bluest galaxies, which have little dust, 
contribute substantially to the 850 /im EBL, despite the 
fact that together they comprise approximately half of 
the sample. 

This result has profound implications for our under¬ 


standing of the star-formation history, lowering previous 
estimates of the high-redshift star formation rate den¬ 
sities by factors of at least two. We analyze the star- 
formation history of our NIR sample using a 20 cm stack¬ 
ing analysis and compare this with the maximum star 
formation rate density at higher redshifts obtained di¬ 
rectly from the submillimeter light. Together these show 
that the total star formation rate density peaks at a red¬ 
shift at or just below one and is roughly fiat at higher 
redshifts. 

The paper is organized as follows. The submillimeter, 
NIR, MIR, optical, radio, and X-ray data are described 
in i|21 The spectroscopic and photometric redshifts are 
discussed in m The use of the NIR, MIR, and radio 
populations to identify the submillimeter background is 
discussed in 0 and the 850 /im EBL identified by the 
NIR-sample is broken down by galaxy flux, color, spec¬ 
tral type, and redshift. The star-formation history is 
described in 0 Our main results are summarized in m 
Throughout the paper, we assume the WMAP cosmol- 
ogy: Hn = 71 km s~^ Mpc“^, Dm = 0.73, and Ha = 0.27 
llBennett et 111 12001 . 

2. THE DATA SAMPLES 

2.1. Submillimeter Data 

The GOODS-N submillimeter map of 
iWang. Cowie. fc Barged ll2004ll is based on jiggle- 
map data taken primarily by our group using the 
SCUBA instrumen t on the JCMT. The data set has also 
been analyzed by IB or vs et al.l ll2003L 120041 . The mo¬ 
saicked submillimeter image has non-uniform 0.4 — 4 mJy 
point-source sensitivity and covers 0.034 deg^ (43% of 
the MIPS area) in the GOODS-N field. Forty-five 
3(7 and s eventeen 4cr submillimeter sou rces are in the 
catalog of fWang. Gowie. fc Barged II2004I1 . We also used 
the submillimeter data to determine the 850 /^m fluxes 
of the various samples. We measured the submillime¬ 
ter fluxes and errors for these samples using optimal 
beam -weighted extractions llWang. Cowie. fc Barged 
120041 throughout the area covered by our submillimeter 
image. 

2.2. Near-Infrared Data 

We carried out deep J and i7-band imaging of the en¬ 
tire GOODS-N region using the Ultra-Low Background 
Camera (ULBCAM) on the Unive rsity of Hawai i 2.2 m 
telescope during 2004 and 2005 llTroiiille et al.ll2?)?)^ . 
ULBCAM consists of four 2kx2k HAWAH-2RG arrays 
llLoose et al ] |200D1 with a total 16' x 16' field of view. 
The images were taken using a 13-point dither pattern 
with ±30" and ±60" dither steps in order to cover the 
chip gaps. The data were flattened using median sky 
flats from each dither pattern. The image distortion was 
corr ected using the ast rometry in the USNO-Bl.O cata¬ 
log llMonet et alJl200,^ . The flattened, sky-subtracted, 
and warped images were combined to form the final mo¬ 
saic with a 20' x 20' area fully covering the GOODS-N 
region. The integration times at each pixel are 19 hours 
in J and 12.5 hours in H, respectively, and the 5a sen¬ 
sitivities are 0.84 /rJy and 2.06 /rJy corresponding to 5a 
AB magnitude limits of 24.1 and 23.1, respectively. A 
more extensive description of the data reduction and a 
detailed analysis may be found in ITrouille et al.l ll2006ll . 
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We ge nerated our source catalo gs with the SExtractor 
package llBertin fc Arnout^IlQQfiH . Because the typical 
seeing was 0"7 and many of the sources appear extended 
in the images, we used auto aperture in SExtractor to 
ensure that the measured fluxes are close to the total 
fluxes. In our stacking analysis, we also consider a cata¬ 
log of all H{AB) < 24 (roughly the Scr limit) sources that 
lie within the submillimeter image. Although this cata¬ 
log will include a small number of false detections 1%, 
inferred from negative sources), this is not a significant 
issue for our stacking analysis, since these sources will 
not affect the signal and will only add a small amount to 
the noise. 

2 . 3 . IRAC Data 

We used the GOODS-N Spitzer Legacy Science Pro¬ 
gram first, interim, and second data relea se products 
(DRl, DR1-I-, DR2; IDickinson et al.lI200fifl . We com¬ 
bined the reduced DRl and DR2 IRAC superdeep im¬ 
ages, weighted by exposure time, to form 3.6, 4.5, 5.8, 
and 8.0 fim images that fully cover the GOODS-N area. 
We again generated the sour ce catalogs in each ban d 
with the SExtractor package llBertin fc Arnoiit^ 1199611 . 
We detected approximately 9800, 8500, 3600, and 3000 
sources at > 5a at 3.6, 4.5, 5.8, and 8.0 /rm, respectively. 
We measured the source fluxes with fixed apertures of 4'.'8 
(3.6 and 4.5 ^m) and 6" (5.8 and 8.0 ^m). These aper¬ 
tures are approximately three times the ~ 177 (3.6 pm) 
to ~ 2" (8.0 pm) EWHM of the point spread function 
(PSF) and are a good compromise between the PSF size 
and the source separation. We applied aperture correc¬ 
tions from the IRAG in-flight PSFs (January 2004) to 
the measured fluxes. The aperture corrections we used 
are consistent with the ones published in the IRAG Data 
Handbook. The corrected IRAG fluxes should be rea¬ 
sonably close to the total fluxes of the sources, because 
the majority of the sources are point-like compared to 
the ~ 2" IRAC PSF. The primary errors in the photom¬ 
etry are caused by the high density of sources, especially 
at 3.6 and 4.5 ^m. In these two bands, the typical dis¬ 
tance between sources is comparable to the PSF, and 
the maps are confusion-limited. Consequently, both the 
background estimate and the aperture photometry are 
highly subject to blending with nearby sources. 

2 . 4 . MIPS Data 

At the longest MIR wavelengths, our sample con¬ 
sists of the MIPS 24 ^m GOODS-N data. We directly 
used the DRl-l- MIPS source list and the version 0.36 
MIPS map provided by the Spitzer Legacy Program. 
The source catalog is flux-limited at 80 qJy and is a 
subset of a more extensive catalog lIGharv et al.ll200(^ . 
With the 0.065 deg^ area coverage, the catalog contains 
1199 24 Atm sources a nd is > 80% complete at 80 ^Jy 
llPapovich et al.ll2004ll . The source positions are based 
on sources detected in the deep IRAC images, and the 
fluxes are derived using PSF fitting. The flux limit and 
accuracy of this catalog are sufficient for our purposes, 
so we did not attempt to generate our own MIPS source 
catalog. 

Nevertheless, because the DR-I- MIPS source list is 
still preliminary, we verified the sources in the list us¬ 
ing the MIPS map. We used a normalized MIPS PSF to 
convolve with the MIPS map and measured the 24 pm 


fluxes at the cataloged positions. These PSF-weighted 
fluxes are mostly consistent with the fluxes in the DR-I- 
source list. However, there are a few sources at the 
edges of the MIPS map, where it is noisier, which have 
low-significance fluxes in our measurements. Seventeen 
of these sources have no obvious counterparts in the 
deep IRAC, NIR, and optical images. These sources are 
likely spurious. Thus, for our stacking analysis, we used 
a restricted area that was fully covered by the Hubble 
Space Telescope ACS GOODS-N observations, since this 
is where nearly all of our submillimeter coverage is. This 
fully avoids the edge problems described above. 

A —0738 offset in declination was applied to t he source 
posit ions to match the radio-frame astrometry llRichardsI 
[ 2 ^. 

2 . 5 . Optical Data 

iCapak et al.1 ll2004ll presented ground-based deep opti¬ 
cal imaging of a very wide-held region around the HDF- 
N. The imaging covers the whole MIPS and IRAC area 
at U, B, V, R, I, z\ and HK' bands. We searched 
fo r counterpart s to th e various samples in the catalog 
of ICaoak et alJ (1200411 using a 1" search radius. This 
search radius closely matches the PSF in the optical and 
the astrometry errors in the optical and MIR. 

Where the images overlapped, we also cross-identihed 
the various samples wi th the ACS GOODS-N catalog of 
iGiavalisco et al.l ll2004ll . 

2.6. Radio and X-ray Data 

We used the 1.4 GHz catalog and image from iRichardsI 
(| 2003 ), which contains sources to a hux limit of 40 /iJy 
(5a), to analyze the radio-selected submillimeter sample. 
We also used the radio data to determine the 20 cm fluxes 
of the various samples. We measured the 20 cm fluxes 
in 3” diameter apertures, adjust ing the normaliz ation to 
match the measured fluxes in the lRichardsI 11200011 catalog 
for the overlapping set of objects. We measured the noise 
level by determining the fluxes at a large number of ran¬ 
dom positions and then measuring the dispersion. The 
la noise is 14 pjy, whic h is almost a fact or of two higher 
than that measured bv iRichardsI ll2000ll . This reflects 
the large aperture used. However, the noise distribution 
measured in this way is well fitted by a Gaussian, with 
an average zero flux level. 

Finally, we use d the Chandra Deep F ield-North (CDF- 
N) 2 Ms catalog llAlexander et al.l2003ll to determine the 
X-ray properties of the various samples and to identify 
sources that contain AGNs. 

2 . 7 . Data Summary 

We summarize the flux limits of the various samples 
in Table ^ We also show the geometry of the various 
data sets schematically in Figure Q] The ground-based 
optical and NIR images cover the entire area shown. We 
mark the region with complete coverage from the ACS 
GOODS-N data with a rectangle and the submillimeter 
region with the contours. The overlap area constitutes 
our core area and covers 106 arcmin^. Essentially all of 
this region lies within 8' (circle) of the X-ray and ra¬ 
dio centers, where the X-ray and radio images have rela¬ 
tively uniform sensitivity. We denote the 4cr submillime¬ 
ter sources with large open squares, the 20 cm sources 
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TABLE 1 

Optical, NIR, and MIR Photometry 


Band 

Sensitivity Limit 

(pJy) 

Telescope 

Reference 

MIPS 24 lira 

80.0 

Spitzer 

a 

IRAC 3.6 /iin 

0.327 

Spitzer 

a 

IRAC 4.5 /iin 

0.411 

Spitzer 

a 

IRAC 5.8 /im 

2.27 

Spitzer 

a 

IRAC 8.0 fim. 

2.15 

Spitzer 

a 

u 

0.052 

KPNO 4 m 

b 

B 

0.063 

Subaru 

b 

V 

0.069 

Subaru 

b 

R 

0.083 

Subaru 

b 

I 

0.209 

Subaru 

b 

z’ 

0.251 

Subaru 

b 

J 

0.839 

UH 2.2 m 

c 

H 

2.06 

UH 2.2 m 

c 


Note. — The MIPS 24 sample is flux-limited and 
complete at 80 fiJy. The median la sensitivity of the MIPS 
map is 6.4 ^Jy. For the rest of the bands, the sensitivity 
limits are 5a limits. References: (a) GOODS Spitzer Le gacy 
Program DRl, DRl-h, and DR2; ('bi lCanak et a1.1120041! : (c) 
this work. 
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Fig. 1. — Schematic layout of the various images. The rectan¬ 
gle shows the deepest portion of the ACS GOODS-N image. The 
contours show the SCUBA image. The circle shows an 8' radius 
around the radio and X-ray centers. The large open squares denote 
the 17 SCUBA sources detected at the 4a level, the solid squares 
denote the 20 cm sources within the 8' radius, and the dots denote 
the 24 fim sample. The total submillimeter area is 125 arcmin^, 
and the area of overlap with the ACS GOODS-N rectangle, which 
we use as our primary area, is 106 arcmin^. The AGS GOODS-N 
rectangle has an area of 144 arcmin^. 


with smaller solid squares, and the 24 /im sources with 
dots. Our core area excludes the more poorly sampled 
regions of the 24 ^m image. 

3. REDSHIFTS 
3.1. Spectroscopic Redshifts 

Intensive spectroscopic redshift surveys have been car¬ 
ried out in the ACS GOODS-N region. We searched 
for spectroscopic redshifts for the sources in I Wirt h et al.l 


(120041 ■ iCowie et al.1 (120041 . and lChaoman et HI (|2005l . 
A substantial number of additional redshifts, which ei¬ 
ther lie outside the ACS GOODS-N region or were ob¬ 
tained from our spectroscopic runs with the Deep Extra- 
galactic Imaging Multi-Object Spectrograph (DEIMOS; 
lEa.ber et 3.11200.*^ on the Keck 10 m telescope subsequent 
to the publication of these papers, were also included. 

3.2. Photometric Redshifts 

We used the U to 8 /im photometry of the sources to 
derive their photometric redshifts. We only considered 
sources that were detected in at least five bands. Gom- 
pared to most optical photometric redshifts, adding the 
deep J and H magnitudes and the Spitzer data has the 
advantage of improving the high-redshift end of the pho¬ 
tometric redshift determinations. The MIR photometry 
probes various spectral features, including the 1.6 /im 
bump caused by the opacity minimum in the stellar at¬ 
mosphere and PAH emission at 5 — 9 prm. In fact, the 
Spitzer photometry alone has been used to derive pho- 
tometric redshifts based on the 1.6 /rm bump (see e.g., 
ISawickilf2002l: lEgami et al.ll2004l . 



Fig. 2.— SED templates used in the training-set photometric 
redshift determinations. Vertical axes are arbitrarily scaled bolo- 
metric fluxes {vFi,). The numerical labels give the SED classes, 
which range from class 1 (elliptical galaxies) to class 7 (very blue 
star-forming galaxies). 

The primary difficulty of incorporating the MIR 
data into the photometric redshift estimation is the 
lack of optical to MIR g alaxy spectrum templates. 
iPerez-Gonzalez et al.l (120051) overcame this problem by 
building “training-set” spectral energy distributions 
(SEDs) from the sources within the MIPS sample it¬ 
self. We also used this method to generate our photo¬ 
metric redshifts. We used just over 1200 galaxies with 
known redshifts and spectral types in the AGS GOODS- 
N region to construct seven templates over the frequency 
range from 6 x 10^^ Hz to 4 x 10^^ Hz. These templates 
range from an elliptical galaxy spectrum to a very blue 
star-forming galaxy spectrum. The seven templates are 
shown in Figure |21 We then made a least-squares fit to 
these templates to determine the photometric redshifts 
and spectral types for the galaxies in each sample. 

The method works extremely well over a wide range of 
redshifts and only fails for a small number of sources. 
In Figure 13 we compare our photometric redshifts 
with the spectroscopic redshifts for the spectroscopically- 
identified H <28 sample in the ACS GOODS-N region. 
There are 1213 sources in this sample with spectroscopic 
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TABLE 2 

Properties of the 850 /xm Sources with 20 cm and 24 Matches 


ID 

^sp 

^ph 

SsSOfiui 

(mJy) 

■S’!.4GHz 

(Mly) 

SsSlim 

(liJy) 

S24pm 

(AiJy) 

Log(5'Hx) 

(ergs cm“^ 

GOODS 850-3 

1.865 

2.71 

7.7 ± 1.0 

151 

14.4 

330 


GOODS 850-6 


2.17 

13.6 ±2.3 

107 

9.5 

185 


GOODS 850-7 

2.578 

3.75 

6.2 ± 1.0 

53.9 

9.8 

313 

-15.10 

GOODS 850-9 

2.490 

2.31 

7.1 ± 1.2 

45.3 

16.0 

235 

-15.04 

GOODS 850-11 


2.67 

10.8 ±2.2 

124 

4.4 

165 

< -15.85 

GOODS 850-15 


2.91 

8.7 ± 2.0 

148 

13.4 

370 


GOODS 850-16 


3.79 

12.4 ± 2.9 

324 

15.4 

267 


GOODS 850-17 

1.013 

1.22 

5.7 ± 1.4 

81.4 

86.7 

724 

-14.65 


redshifts that are not saturated in the optical {z' > 19), 
and 1134 of these have statistically acceptable fits to the 
templates at some redshift from z = 0 to 4. The remain¬ 
ing sources either have unusual SEDs or are blended with 
a neighbor. For the 1134 sources with both spectroscopic 
and photometric redshifts shown in Figure 13 there are 
only a couple of seriously discrepant sources, and while 
the scatter becomes larger at the higher redshifts, the 
method robustly places nearly all the sources in the cor¬ 
rect redshift range. 



Spectroscopic Redshift 


Fig. 3.— Comparison of our photometric redshifts, as deter¬ 
mined from the training-set templates, with the spectroscopic red¬ 
shifts for the spectroscopically-identified FT-band sample in the 
ACS GOODS-N region. 


4. IDENTIFICATION OF THE 850 MICRON EBL 

4.1. Direct Identification of Bright Suhmillimeter 
Sources 

Many of the bright (> 2 mJy) submillimeter sources 
that are directly detected in the 850 pm blank-field 
SCUBA images can be localized by their radio emis¬ 
sion, as discussed in the introduction. However, since 
the 24 pm-band traces the short wavelength end of the 
FIR emission, the 24 pm sample may provi de additional 
count erparts to the brig ht SCUBA sources llEgami et al.l 
l2nnl llvisonTTaniWil . 

Our 850 pm SCUBA map covers 0.034 deg^, and there 
are 17 4g bright SCUBA source s in th is area in the cat¬ 
alog of lWang. Cowie. fc Bargeil ll2004H . whose notation 


we shall follow. Three of the bright SCUBA sources 
(GOODS 850-2, 850-8, and 850-14) have no 24 pm 
counterparts, even within a very wide 8" search ra- 
dius. These sources also do no t have radio counterparts 
ijWang. Cowie. <C Ra,rgerl|2flfl4jl . The remaining 14 have 
a total of twelve 24 pm sources within a 6" search radius 
and eighteen 24 pm sources within an 8" search radius. 
Given the number density of 24 pm sources, we expect 
~ 3 false matches at 6" and ^ 5 false matches at 8", 
so we are identifying some 9 — 13 of the remaining 14 
sources, or about 50 — 75% of the total sample. 

We can compare this with the counterpart 
identifications made using the 20 cm data by 
IWang. Cowie. fc Barged ll2004H . Eight of the 17 do- 
bright SCUBA sources have 20 cm counterparts within 
6", and we expect all of these to be real, given the surface 
density of radio sources. Ten have counterparts within 
8", including one source with two radio counterparts; 
we expect one false match at this radius. These results 
suggest that the 24 pm data may be slightly better at 
picking up the bright SCUBA sources, but the difference 
is not statistically significant. All eight radio sources 
with bright SCUBA counterparts within 6" are seen at 
24 pm, and the 24 pm data then picks out a further four 
sources, of which three are expected, on average, to be 
false. 

The redshifts (spectroscopic and photometric) and 
properties of our eight securely identified bright SCUBA 
sources are summarized in Table 13 The whole 24 pm 
sample has a median 24 pm flux of 152 pjy. The SCUBA 
sources in Table |3 have a median 24 pm flux of 300 pJy 
and are weighted to the high end of the 24 pm fluxes. 
They are mostly at high redshifts {z = 2 — 3), with a 
median redsh ift of z = 2.5 , and are consistent with be¬ 
ing ULIRGs. Ilvison et alJ ll2002ll presented a sample of 
30 radio-identified SCUBA sources with 850 pm fluxes 
> 8 mJy. Using the millimetric redshift technique, the-y 
found a median redshift of z = 2.4. iChaoman et alJ 
Il2nn3all presented a spectroscopic redshift sample of 
10 radio-identified bright (> 2 mJy) SCUBA sources, 
and IChaoman et al.l ll2005fl expanded that sample to 73 
sources. Their typical redshift range isz = 1.7 — 2.8, 
with a median redshift of z = 2.2. The redshift distribu¬ 
tion of our radio-identified bright SCUBA sample is fully 
consistent with the results of the above groups. 

4.2. A Stacking Analysis of the Suhmillimeter EBL 
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TABLE 3 

850 /xM Stacking Analyses 


Waveband 

FlUX77T,i7T, 

(/^Jy) 

N 

Area 

(arcmin^) 

< Ssso > 
(mJy) 

I. 

(Jy deg-2) 

clean 

(Jy deg-2) 

ACS 0.45 fim 

0.05 

6868 

106 

0.10 ±0.017 

24.8(15.2) ±4.1(6.0) 

14.3(19.3) ±4.0(6.0) 

ACS 0.8 fim 

0.1 

7826 

106 

0.10 ±0.016 

26.8(16.3) ±4.3(6.4) 

15 . 4 ( 20 . 9 ) ±4.2(6.3) 

ULB 1.6 fim 

1 

3094 

122 

0.20 ±0.03 

18.8(21.1) ±2.3(3.5) 

10 . 4 ( 17 . 7 ) ±2.3(3.5) 

IRAC 3.6 fim 

0.3 

5245 

106 

0.11 ±0.0066 

19.6(16.7) ± 3 . 4 ( 5 .!) 

8.2(15.6) ±3.4(5.1) 

IRAC 8.0 fim 

2.0 

1587 

106 

0.33 ±0.03 

18.0(l9.0) ± 1. 9 ( 2 . 8 ) 

9 . 1 ( 13 . 8 ) ± 1 . 9 ( 2 . 8 ) 

MIPS 24 pm 

80 

493 

106 

0.66 ±0.06 

11.4(9.1) ± 1.1(1.5) 

5.9(5.6) ± 1.1(1.6) 

VLA 20 cm 

40 

101 

122 

1.31 ±0.13 

4.0(5.1) ± 0 . 40 ( 0 . 6 ) 

1.5(2.7) ±0.4(0.6) 

JNIR 1.6, 3.6 fim 

1.8 

3121 

106 

0.20 ±0.025 

20.7(20.9) ±2.6(4.6) 

11.4(16.1) ±2.6(4.0) 


Note. — The N values are the numbers of sources with errors less than 4 mjy that were used in 
the stacking analyses. 


We cannot directly analyze the submillimeter sources 
that are too faint to be detected in blank-field SCUBA 
surveys, but we can use our galaxy samples at other 
wavelengths to study these sources using a stacking tech¬ 
nique. This technique was first used bv IPeacock et al.1 
(l2000t) to study the submillimeter properties of Lyman- 
break galaxies. Here we use our SCUBA map to deter¬ 
mine the average 850 fim properties of the sources in 
a given sample, and hence the amount of 850 fim light 
they produce. Although such an analysis provides very 
little information on the individual source properties, it 
can show which of the galaxy populations give rise to the 
submillimeter light. Since we can also break the samples 
down by galaxy flux, color, spectral type, or redshift, we 
can determine the properties and redshift distribution of 
the class(es) of sources that is (are) producing the bulk 
of the 850 ^m EBL. 

One difficulty with trying to assess what fraction of 
the 850 p,m EBL a class of sources produces is that 
the absolute normalization of the 850 /im EBL is some¬ 
what uncertai n. The observed 85 0 /rm EBL is given as 
31 Jy deg~^ i niPuget et ahl lIlQQfifl and as 44 Jy deg“^ in 
iFixsen et alJ l|199S|l . The substantial difference between 
the two estimates is a consequence of different correc¬ 
tions for foreground emission, and it is unclear which is 
the better estimate. We therefore compare with the full 
range in our subsequent analysis. 

The possible source populations that are contributing 
the bulk of the 850 /im EBL are also constrained by 
the deepest sub millimeter number counts . The 850 /tm 
number counts o flCowie. Barger, fc KneibI ll2002ll showed 
that the number of submillimeter sources with fluxes 
> 0.5 mJy is approximately 2.5 x 10'^ deg“^. This is 
therefore a rough minimum on the density of any sam¬ 
ple that seeks to explain the 850 /im EBL. Simply based 
on such number densities, it is clear that there are not 
enough 20 cm selected sources to account for the light, 
and there are only marginally enough 24 /tm sources (ap¬ 
proximately 1.8 X 10'^ deg“^ above 80 /tJy). In other 
words, the current 20 cm and 24 /tm samples are not 
deep enough to have reached the fainter submillimeter 
sources, while the NIR and MIR samples are deep enough 
to have surface densities that are substantially above the 
required value. 

We first computed the contributions to the 850 /tm 
EBL from seven samples: the iL-band, 3.6 /tm, 8 /tm. 


24 /tm, and 20 cm selected samples, as well as the B 
and /-ban d selected sample s from the ACS GOODS-N 
catalogs of iGiavalisco et aO ll2004ll . Eor each source in 
each sample, we determined the beam-we ighte d 850 /tm 
flux and noise from our SCUBA map (see EU- We then 
measured the error-weighted average 850 /tm flux for all 
of the sources in the sample that had errors less than 
the chosen cut value of 4 mJy. We also determined the 
area where each sample had submillimeter measurements 
of this sensitivity. The EBL contribution is then given 
by the product of the number of sources and their error- 
weighted mean divided by the observed area. The results 
are almost identical for cut values other than 4 mJy. The 
contribution of each population to the 850 /rm EBL is 
summarized in Tabled where we give the sample wave¬ 
length and limiting flux, the area covered, the number 
of sources, the error-weighted mean 850 /tm flux for each 
source, the 850 /tm EBL contribution with its la er¬ 
ror, and the 850 /tm EBL contribution measured from 
a map where 4tT 850 /tm sources are CLEANed (see be¬ 
low). Since the error weighting may result in a very small 
region dominating the signal, we also computed the un¬ 
weighted 850 /tm signal over the same area. This is given 
in parentheses in the final two columns of the table. 

Our results confi rm the MIR detections made by 
iSerieant et al.1 l|2004ll with the Spitzer Early Release Ob¬ 
servations and the SCUBA 8 mJy survey, but at substan¬ 
tially higher significance levels. As an example, there are 
493 24 /tm sources with 850 pm errors less than 4 mJy 
in the SCUBA map. The mean 850 /tm flux of these 
sources is 0.66 ±0.06 mJy. This is a substantial improve¬ 
ment in s ignal to noise over th e value of 0.30 ± 0.24 mJy 
found bv ISerieant et al.1 l|2004ll due to the large number 
of 24 /tm sources in the GOODS-N field. We obtain a 
consistent average flux of 0.53 ± 0.08 mJy from the un¬ 
weighted sample, which shows that the lower error re¬ 
gions are not dominating the signal. The total stacking 
contribution of the 24 pm sources to the 850 pm EBL 
is 11.4 ± 1.1 Jy deg“^, or about 25 — 35% of the total 
850 pm EBL. 

As expected, given that they do not have sufficiently 
high number densities, the 20 cm and 24 pm samples 
only identify a fraction of the 850 pm EBL. (The radio 
sample identifies 4.0 ± 0.40 Jy deg“^, or about 10 — 13% 
of the 850 pm EBL.) It is clear that much deeper samples 
would be required at both of these wavelengths to sub- 
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stantially identify the 850 ^m EBL. To try to see how 
much fainter we would have to go at 24 ^m, we gen¬ 
erated a fainter 24 ^m source catalog using SExtractor 
and excluded the original 1173 MIPS sources. There are 
731 3(7 MIPS sources in this catalog, and 338 of them 
are in the SCUBA map. Their median 24 /rm flux is 
46 /rJy, and the faintest ones are ~ 10 /rJy. These faint 
24 ^m sources possess a weak 850 ^m stacking signal of 
1.73 ± 1.32 Jy deg“^. We note that the completeness of 
the 24 /rm catalog decreas es rapidly from ~ 1.0 at 80 /rJy 
to ~ 0.1 at 40 /rJy te.g.. iPapovich et al.ll2004ll . There¬ 
fore, the faint 24 /im sources in the 10 — 80 /iJy range 
might make a significant contribution (> 10 Jy deg“^) 
to the 850 pm EBL. However, because of the low statis¬ 
tical significance of the stacking signal, we cannot show 
that this is the case with the current data. To provide 
a definitive answer, we would need a substantially com¬ 
plete 24 pm catalog to about a 10 pJy flux. 

By contrast, the 8 pm sample provides a much more 
substantial identification of the 850 pm EBL (18.0 ± 
1.9 Jy deg“^), as do all of the shorter wavelength bands. 
In fact, at first sight, the optical wavelength bands give 
the most complete identification of all. However, on 
closer inspection, much of this signal comes from the very 
small HDF-N region, where the submillimeter errors are 
very low. Here the slightly smaller EBL determined from 
the unweighted average may be a more representative 
value. 

The signal from each of the samples is, of course, heav¬ 
ily overlapped, and we next compared the samples to see 
how much additional signal each added. To do this, we 
took each sample in turn as our primary sample and then 
measured signals from the residual sources in each of the 
other samples after each of the primary sample’s over¬ 
lapping sources were excluded. That is, we measured the 
signal from the other samples, after excluding all sources 
which were already present in the primary sample. We 
were able to determine which of the samples were most 
complementary based on this procedure. As a result of 
our analysis, we formed a sample from the combination 
of all sources with iJ-band or 3.6 pm fluxes greater than 
1.8 pJy. We restricted the sample to the ACS GOODS-N 
region. We refer to this sample as our NIR sample, and 
we summarize its properties in the final line of Table El 
Physically, since the two bands bracket the peak flux in 
most galaxy SEDs over roughly the z = 0 — 3 range, this 
sample is choosing nearly all of the galaxies with peak 
observed fluxes above the 1.8 pJy cut. 

Erom our NIR sample, we find a contribution of 20.7 ± 
2.6 Jy deg“^, or 50 — 70% of the 850 pm EBL. Perhaps 
even more importantly, our NIR sample appears to con¬ 
tain nearly all of the EBL that can be measured from the 
remaining samples. In other words, when our NIR sam¬ 
ple is excluded, none of the residual 8 pm, 24 pm, or ra¬ 
dio sample sources gives a signal greater than 1 Jy deg“^. 
The residual B and /-band sample sources give a signal 
of a few Jy deg“^ in the weighted samples but a null sig¬ 
nal in the unweighted samples. We therefore adopt our 
NIR sample as our primary sample for further analysis. 

There are 3121 sources in our NIR sample. Of these, 
2415 are in the original // < 24 3tT catalog, and a further 
374 have //-band magnitudes brighter than 24, reflect¬ 
ing the incompleteness of this catalog in the H = 23 — 24 
range. The remaining 332 sources in our NIR sample 


would not have been selected in a complete // = 24 sam¬ 
ple and have only been detected using the 3.6 pm cata¬ 
log; they give a signal of 4.8 ± 0.9 Jy deg“^. Thus, the 
combined catalog gives an improved identification of the 
EBL relative to either catalog alone. 

We visually inspected the 3.6 pm image for all of the 
sources without optical or NIR counterparts to check 
that these were not spurious. About 30 sources are 
suspect, either because they are not clearly seen in the 
3.6 pm image, or because they are contaminated by a 
neighbor. However, excluding these sources has no effect 
on the measured signal, and we conclude that there is a 
significant contribution from the 3.6 pm sources in the 
sample. 

To check the robustness of our results, we performed a 
number of tests. First, we summed the total SCUBA im¬ 
age to find the total flux of our SCUBA map and found 
the result to be strictly zero. This is as expected, since 
each positive source in the SCUBA map has two negative 
50% sidelobes. The zero sum of our SCUBA map indi¬ 
cates that there was good sky subtraction during the data 
reduction. Because of the zero sum of the submillimeter 
map, any random populations will have zero stacking 
fluxes, and a positive stacking signal will indicate a cor¬ 
relation between that population and the submillimeter 
sources. To test this and the assigned statistical errors, 
we measured fluxes for large numbers of random posi¬ 
tions in the SCUBA image and analyzed these in the 
same fashion as the real samples. The results were fully 
consistent with expectations in both the average signal 
and the statistical spread. 

The third test was a Monte Carlo simulation. We 
used the 24 pm source catalog and randomized its as¬ 
trometry to create simulated 850 pm sources. The 
850 pm fluxes of the simulated sources were derived 
from conversions based on M82 and Arp 220. The sim- 
ul ated sources were added into the “true noise” map 
of rWang. Cowie. fc BargeJ ll2004H . which has both the 
bright SCUBA sources and the faint confusion sources 
removed, such that only noise is left. When we mea¬ 
sured the 850 pm fluxes from the simulated maps, we 
randomly offset the measured positions from the source 
positions with up to 2" rms position errors, enough to 
account for the pointing errors of the submillimeter tele¬ 
scope and the astrometry errors in the 24 pm catalog. 
The average stacking flux was measured in 100 such re¬ 
alizations using the same methodology as was used in 
analyzing the true submillimeter image and fully recov¬ 
ered the known input fluxes. In summary, our Monte 
Carlo simulations and the zero sum of the SCUBA map 
show that the stacking flux is an unbiased estimate of the 
850 pm EBL and that the assigned errors are realistic. 

As our fourth test, we split the NIR sample into two 
parts: one corresponding to sources where the measured 
submillimeter errors lay between 0.1 and 1.5 mJy, and 
the other to sources where the errors lay between 1.5 
and 4 mJy. The first sample contained 826 sources and 
gave an 850 pm EBL contribution of 20.1 ±3.2 Jy deg“^. 
The second sample contained 2295 sources and gave an 
850 pm EBL contribution of 21.8 ± 4.5 Jy deg“^, show¬ 
ing that an equivalent signal can be obtained from two 
samples with very different sensitivities. 

Correlations in the target sample can result in an up¬ 
ward bias to the signal if the submillimeter flux comes 
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from overdense regions, since, in this case, we over- 
count the mean flu x associated with each source (e.g., 
iSerieant et al.ll2?)?)^ . Given the expected angular corre¬ 
lations in the samples, we do not expect this effect to be 
large, but it may be present at some level. We tested this 
by measuring the average number of galaxies lying within 
7” of a given galaxy. Since this is the half-width of the 
850 fiia beam, we expect that only galaxies within this 
area could be substantially contaminated. For the NIR 
sample described above, which has a surface density of 
1.1 X 10^ deg-2, we found an average number of galaxies 
of 1.1 within this radius. This is actually smaller than 
the expected number for a random distribution (1.18), 
so there is very little correlated signal, and the correla¬ 
tion effects on the stacking signal are small. Since we 
shall subsequently restrict the sample even further, we 
ran this same test on that restricted sample and found a 
similar result. 

Other groups te.g-.. ISerieant et al.l2004ll perform stack¬ 
ing analyses on maps where known submillimeter sources 
are removed. As our final test, we also removed all known 
submillimeter sources detected at or above the 4 tr level 
from the map and reran our stacking analyses. To do 
this we subtracted (CLEANed) the fitted PSFs of all of 
these sources from the image. Because our noise esti¬ 
mate does not include the confusion noise (e.g., Cowie et 
al. 2002; Wang et al. 2004), our 4 a cut is comparable 
to the 3.5 tr cut used by most other groups. We note 
that such a sigma cut has different flux limits in differ¬ 
ent areas of our map. The totol flux removed from the 
map is only 5.6 Jy deg“^ but the sources removed from 
the deepest areas indeed correspond to a total surface 
brightness of ~ 10 Jy deg“^ according to the number 
counts in Wang et al. (2004). The removal therefore af¬ 
fects the weighted contribution to the EBL more than 
the unweighted contribition and the later value is com¬ 
parable to the reduction in the weighted signal while the 
former matches the change in the unweighted signal. As 
a consequence of this complex biasing effect from the 
non-uniform sensitivities, we decided to use only our un- 
CLEANed results. However, for completeness, we also 
include our CLEANed results in Table 3. 

4.3. The 850 fim EBL versus Galaxy Flux, Color, and 
Spectral Type 

We can further subdivide the contributions to the 
850 /rm EBL by galaxy flux, color, and spectral type 
to determine the properties of the sources giving rise to 
the light. 

We show the contributions to the 850 pm EBL from 
the iL-band sample versus the iL-band flux in Figure 01 
and from the 3.6 pm sample versus the 3.6 pm flux in 
Figure 01 In each case, we denote with solid squares the 
contributions in half dex flux intervals, starting at the 
limiting flux of the sample. The total contribution of 
the entire sample is given in the upper right corner. The 
individual source fluxes in mJy are denoted by dots; 4cr 
measurements are denoted by plus signs. For the lowest 
flux interval in the i7-band sample, we also show the EBL 
value corrected for the incompleteness in the iL-band cat¬ 
alog (open square). For both samples, the contributions 
are dropping at the faintest end, which would be con¬ 
sistent with the onset of convergence to the asymptotic 
value. The peak contribution comes near a flux of 6 pJy 


in the iJ-band and near a flux of 18 pJy in the 3.6 ^m- 
band. Expressed in magnitudes, the peak contribution 
to the measured signal is coming from sources with NIR 
magnitudes of 21 — 22 (AB). 



Flux (microJy) 


Fig. 4.— Contributions to the 850 fim EBL from the Lf-band 
sample vs. the iL-band flux. The solid squares show the contri¬ 
butions from each half dex flux interval with Icr error bars. The 
number in the upper right corner gives the total contribution. The 
open square shows the contribution from the lowest flux interval 
corrected for the incompleteness in the iL-band catalog. The dots 
denote the measured fluxes of the individual sources for a y-axis in 
mJy units, and the plus signs denote sources with 4cr measurements 
at 850 pm. 



Fig. 5.— Contributions to the 850 pm EBL from the 3.6 pm sam¬ 
ple vs. the 3.6 pm flux. The solid squares show the contributions 
from each half dex flux interval with Icr error bars. The number 
in the upper right corner gives the total contribution. The dots 
denote the measured fluxes of the individual sources for a y-axis in 
mJy units, and the plus signs denote sources with 4cr measurements 
at 850 pm. 

We show the contributions to the 850 pm EBL from 
the iL-band sample versus the I — H color in Figure 01 
and from the 3.6 ^m sample versus the H — 3.6 pm 
color in Figure [ 7 | In both cases, the light comes from 
the redder sources in the sample. In the iJ-band sam¬ 
ple, the I — H color weighted by the submillimeter light 
contribution is 1.5, while the mean color of the whole 
sample is 1.0. This arises because there is a higher 
mean submillimeter signal per object (open diamonds) 
in the red-colored sources. Similarly, the 3.6 ^m sam¬ 
ple has a submillimeter-weighted color — 3.6 pm of 
1.2, as compared to a mean color of the whole sample 
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I-H (AB) 


Fig. 6 .— Contributions to the 850 //m EBL from the ff-band 
sample vs. the I — H color in AB magnitudes. The solid squares 
show the contributions from each color interval with Icr error bars. 
The open diamonds show the mean 850 pm fluxes of the sources 
for a y-axis in mjy units. 



H-3.6micron (AB) 


Fig. 7.— Contributions to the 850 pm EBL from the 3.6 pm 
sample vs. the H — 3.6 pm color in AB magnitudes. The solid 
squares show the contributions from each color interval with Icr 
error bars. The open diamonds show the mean 850 pm fluxes of 
the sources for a y-axis in mJy units. 


of —0.3. The cumulative signal in the i?-band sample is 
3.4 ± 1.1 Jy deg“^ ior I — H > 2 and 14.5 ± 2.1 Jy deg“^ 
ior I - H > 1. 

This biased contribution of red galaxies to the submil¬ 
limeter light is a known result from stacking analyses car¬ 
ried out using submillimeter measurements of extremely 
red o bjects in both lensed cluster and blank-field sur¬ 
veys iWe hner. B arg-er. fc Kneib ll20f)^ iWebb et al.ll2004 
iKnudsen et alJ l2005D . 

However, the observed-frame colors are a function of 
galaxy type, reddening, and redshift, and we may make 
a much more powerful analysis using the photometric 
redshift determinations, which separately yield the red- 
shift and the rest-frame SED for each source. In Fig¬ 
ure |H1 we show the 850 ^m EBL contribution from the 
galaxies in the NIR sample divided into the seven SED 
classes used in our training-set photometric redshift anal¬ 
ysis (see Fig. El. We assign unidentified sources an SED 
class of zero, and we exclude spectroscopically-identified 
stars. The vertical lines split the sample into four bins: 
unidentified sources (class 0), elliptical galaxies (class 1), 
intermediate spiral galaxies (classes 2 — 5), and very blue 
star-forming galaxies (classes 6 and 7). We have printed 
directly on the figure the number of sources in each of 


the four bins. 



SED Class 

Fig. 8 .— Contributions to the 850 pm EBL vs. the SED class. 
Unidentified sources are placed in class 0, and spectroscopically- 
identified stars are excluded. The solid squares show the contribu¬ 
tions from each class with Icr error bars. The vertical lines divide 
the sources into 4 bins, with the number of sources in each bin 
printed at the top of that bin. 
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SED Class 

Fig. 9.— Contributions to the 20 cm EBL vs. the SED class. 
Unidentified sources are placed in class 0, and spectroscopically- 
identified stars are excluded. The solid squares show the contribu¬ 
tions from each class with Icr error bars. (These are smaller than 
the symbol size.) The vertical lines divide the sources into 4 bins, 
with the number of sources in each bin printed at the top of that 
bin. 



Neither the elliptical galaxies nor the very blue star¬ 
forming galaxies give a significant signal. The 1213 
galaxies in these three classes have an error-weighted 
850 ^m flux of —0.03 ± 0.04 Jy, presumably reflecting 
the lack of star formation in the elliptical galaxies and 
the absence of dust in the blue galaxies. Nearly all of 
the 850 fj^ia signal comes from the intermediate spiral 
galaxies, which have an error-weighted 850 /rm flux of 
0.40 ± 0.03 Jy. There is a small contribution to the 
850 fim EBL from the unidentified sources, which have 
an error-weighted 850 fim flux of 0.37 ± 0.14 Jy. 

The same selection appears if we consider the 20 cm 
or 24 fim properties of the NIR sample as a function of 
SED class. In order to determine the average 20 cm flux 
of the NIR sample, we performed a stac king analysis us¬ 
ing the 20 cm image of lRichar'd^ 11200011 . (For the radio 
image, since the errors are constant, the average signal 
is the same as the error-weighted signal.) As with the 
submillimeter analysis, we measured the 20 cm fluxes for 
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SED Class 

Fig. 10.— Fraction of the sources detected above a 24 /im flux of 
80 fiJy vs. the SED class. Unidentified sources are placed in class 
0, and spectroscopically-identified stars are excluded. The vertical 
lines divide the sources into 4 bins, with the number of sources in 
each bin printed at the top of that bin. 


each of the sources in the NIR sample (see the 

average 20 cm flux per source is 5.7 ± 0.34 ^Jy. In com¬ 
puting the 20 cm EBL, we eliminated the small number 
of 20 cm sources brighter than 300 /xJy in the held, since 
none of these are directly detected submillimeter sources. 
The 20 cm EBL is dominated by the fainter sources, and 
this cutoff makes only a small change in the results. In 
Figure El we show the 20 cm EBL contribution from the 
galaxies in the NIR sample divided into the same seven 
SED classes. Here, too, nearly all of the radio light comes 
from the intermediate SED classes (2 — 5). These inter¬ 
mediate spiral galaxies have a signal of 9.1 ± 0.25 /xJy. 

In Figure nni we show the fraction of sources in each 
SED class that are detected above 80 fiJy at 24 ^m based 
on the 24 /im catalog. Nearly all of the 24 /rm sources 
with fluxes above 80 /rJy also lie in the intermediate SED 
classes. 

Given that the elliptical galaxies and the very blue 
star-forming galaxies are only a source of noise, we 
now remove these objects and the spectroscopically- 
identified stars from the NIR sample, leaving us with 
1770 sources. We refer to this sample as our core NIR 
sample. The 850 /rm EBL for these sources is plotted 
versus the greater of the i7-band or 3.6 /rm flux in Fig¬ 
ure m The total contribution from this final sample is 
24.0±2.0 Jy deg“^, or about 54—77% of the total 850 /rm 
EBL. 

4.4. The Redshift Distribution of the 850 frm EBL 

We may now use the photometric and spectroscopic 
redshifts to determine where in redshift space the 850 /rm 
EBL arises. Nearly half of the sources (1478) in the 
full NIR sample have spectroscopic redshifts. Combin¬ 
ing these with the photometric redshifts increases the 
identification to 3020. Only 101 sources are blended, too 
faint, or too peculiar in their SEDs to be identified. 

In Figure [El we show the 850 p,m EBL that arises 
in the core NIR sample (that is, the intermediate SED 
classes plus the unidentified sources only) divided by red¬ 
shift interval. Here the solid squares denote the con¬ 
tributions from the sample with either photometric or 
spectroscopic redshifts, and the open squares denote the 
contributions from the sample with spectroscopic red- 
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Fig. 11.— Contributions to the 850 /rm EBL from the unidenti¬ 
fied (class 0) and intermediate spectral type (classes 2 — 5) sources 
(i.e., our core NIR sample) vs. the greater of the H or 3.6 pm band 
fluxes. The solid squares show the contributions from half dex flux 
intervals with Itr error bars. The number in the upper right corner 
gives the total contribution. The dots denote the measured fluxes 
of the individual sources for a y-axis in mjy units, and the crosses 
denote sources with 4(7 measurements at 850 pm. 


shifts only. Open diamonds show the average flux per 
source with a y-axis in mJy. The EBL from the core 
NIR sample is dominated by low-redshift sources. In¬ 
deed, 14.0 ± 1.6 Jy deg“^ of the EBL comes from below 
z = 1.5, implying that about half of the 850 ^m EBL 
is originating at these low redshifts. This is in striking 
contrast to the redshift distribution of the bright SCUBA 
sources seen at higher fluxes and identihed using their ra¬ 
dio counterparts. It appears that at lower submillimeter 
fluxes, there is a substantial contribution to the 850 /im 
EBL from galaxies at much lower redshifts than is the 
case at the higher submillimeter fluxes. 



Fig. 12.— Contributions to the 850 /rm EBL from our core NIR 
sample with spectroscopic or photometric redshifts vs. redshift. 
The solid squares show the contributions from each redshift interval 
with 1(T error bars. The open squares show the contributions if 
we only consider sources with spectroscopic redshifts. The open 
diamonds show the mean fluxes of the sources for a y-axis in mJy 
units. 


5. THE STAR-FORMATION HISTORY 

5.1. The Star Formation Rate Density from the Core 
NIR Sample 


FLUX (mJy) 
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Subject to the assumed stellar mass function, the star- 
formation rates (M in Mq yr“^) of sources can be esti¬ 
mated from their infrared luminosities. In our analysis, 
we use the formula M = 1.7 x l|Kennicuttl 

[T99|). This is very similar to the value of M = 1.5 x 
derived by iBarger. Cowie. fc Richard^ 

lf200nll . 

The FIR luminosity can, in principle, be estimated 
from the 3.6 — 24 /rm flux, and this has been done to 
estimate the star-formation history of the 24 fim popu¬ 
lation. However, the conversion from 3.6 — 24 /im flux to 
FIR luminosity is complex, so this method is relatively 
uncertain, and we do not follow it here. 



Fig. 13.— Contributions to the 20 cm EBL from our core NIR 
sample with spectroscopic or photometric redshifts vs. redshift. 
The solid squares show the contributions from each redshift interval 
with 1(T error bars. 

The direct conversion of an 850 /im flux to a FIR lu¬ 
minosity is probably relatively robust for high-redshift 
(z > 1) ULIRGs, and we use this in ^5.21 to estimate the 
star formation rate density (SFRD) for those sources. 
However, as we saw in E3I much of the submillimeter 
EBL that we have been able to identify from the core 
NIR sample is at lower redshifts, and the submillimeter 
flux to FIR luminosity conversion for these lower lumi¬ 
nosity sources may have a much wider range. Thus, in 
order to compute the star-formation history of the core 
NIR sample, we use the 20 cm fluxes with their robust 
conversion of radio power to total luminosity. 

In Figure [CT we show the 20 cm EBL versus redshift. 
We see a strong peak at a redshift of just below one, 
which then trails down to near zero beyond z = 2.6. 
We also generated the same plot for a large number of 
samples of randomized positions in the field containing 
the same number of sources. The random realizations 
average to zero, are fully consistent with the range ex¬ 
pected from the statistical noise, and show no systematic 
effects. In addition, we computed the 20 cm EBL with 
redshift for only the sources in the iH.ichardsI l|200nll cata¬ 
log with its limiting flux of 40 /iJy. As might be expected, 
since this sample is considerably brighter in flux than the 
stacking sample, we found a lower peak redshift, with a 
value of z 0.5. 

We note that Figures El and na imply an unusually 
large submillimeter-to-radio flux ratio at z ~ 1. This is 
consiste nt with cool dust in local low-lu minosity galax¬ 
ies (e.g., IVlahakis. Dunne, fc Ealesll2005ll or with a great 
amount of high-redshift background light being lensed by 


low-redshift sources le.g.. lAlmaini et al ] l200l . 

We next translated the 20 cm EBL into a bolometric 
luminosity density for each redshift interval. We used 
the center redshift of the interval to compute the radio 
power per unit area, which we then converted t o a total 
lumin osity using the FIR-radio correlation (e.g., iCondorJ 
IT^ . We next computed the volume in each redshift bin 
for the unit area, thereby forming the total luminosity 
density per unit volume at that redshift. 

In Figure El we show the SFRD of the core NIR sam¬ 
ple obtained by converting the lu minosity density to a 
SFRD using the iKenniciii^ (j1 relation. The SFRD 
(solid squares) shows a rapid rise to z = 0.8, and then it 
flattens at higher redshifts. Since some part of the light 
may arise from AGNs rather than star formation, we also 
show in Figure m the SFRD from only the sources that 
have X-ray luminosities less than 10"^^ ergs s“^ (open 
diamonds). This has only a small effect, so unless the 
correction for X-ray obscured AGNs is much larger, we 
can ignore the AGN contamination at our current level 
of accuracy. 

At each redshift, the measured SFRD is a lower bound 
on the total submillimeter SFRD, since the core NIR 
sample is not a complete mapping of the submillimeter 
star-formation history. In particular, it may be biased to 
lower-redshift sources and preferentially miss the higher- 
redshift star formers. In asa we assess the maximum 
corrections that are possible from the residual 850 /im 
light that the core NIR sample omits. 



Fig. 14.— Submillimeter SFRD vs. redshift. The solid squares 
show the SFRD derived from the 20 cm EBL contributions of our 
core NIR sample. (The 1 a errors are smaller than the symbol 
sizes.) The open diamonds show the same results when sources 
with a 2 — 8 keV luminosity above 10"^^ ergs s~^ (sources con¬ 
taining AGNs) are excluded. The solid triangles show the SFRD 
computed using the 850 /im EBL. The dashed horizontal line shows 
the SFRD computed assuming all of the sources with 850 /im fluxes 
above 4 mjy lie at 2 = 1 — 3, as suggested by the bright source 
identifications. The rectangular region denotes the SFRD from the 
remaining submillimeter EBL that is not accounted for by our NIR 
sample, assuming that it also lies in the redshift interval 2 = 1 — 3. 
The range corresponds to the uncertainty in the 850 /im EBL. The 
maximum total submillimeter SFRD in this redshift range is then 
the sum of the rectangle and the measured points. 


5.2. Limits on the Star Formation Rate Density at 
z= 1-3 
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For the higher-redshift sources, where the galaxies are 
near-ULIRGs, we may make a direct estimate of the 
SFRD from the submillimeter light. We again assume 
that the luminosities are dominated by star formation. 
The infrared luminosities of the submillimeter sources 
can be estimated from their 850 ^m fluxes, redshifts, a 
plausible dust temperature (Td), and a dust emissivity 
index (/3) in the submillimeter. To make this conver sion, 
we adopted the dust model in lYiin fc Carillil II2002I1 and 
their values of = 58 ± 9 and f3 = 1.32 ± 0.17 based 
on luminous and ultraluminous starbursts. This conver¬ 
sion is very similar to the one obtained for Td = 47 and 
[3 = 1.0 (Arp 220), which has often bee n used in the past 
('e.g.- lBarger. Cowie. fc Richardsll20(inll . 

We first computed the SFRD from the 850 /rm EBL 
determined as a function of redshift for the core NIR 
sample. We computed the FIR luminosity per unit area 
in each redshift bin using the adopted dust model. Next, 
we computed the volume corresponding to the redshift 
interval. Finally, we converted t he FIR luminosit y per 
unit volume to a SFRD using the iKennicuttI lll998ll rela¬ 
tion. The result is plotted in Figure d as the solid tri¬ 
angles. At z > 1, these results agree strikingly well with 
the 20 cm-determined SFRD shown by the solid squares. 
This is not a foregone conclusion, since we are using two 
independent data sets for the computation and two dif¬ 
ferent (though hopefully consistent) calibrations for the 
conversion to luminosity. (The FIR-radio correlation in 
the ca se of the 20 cm conversion, and the lYiin fc CarillH 
(l2003l dust model for the submillimeter conversion.) At 
z < 1, the 850 /im-determined SFRD is higher than the 
20 cm determined SFRD, which is expected, since the 
adopted luminosity conversion may no longer be appro¬ 
priate if the sources are no longer near-ULIRGs. 

We next computed the SFRD assuming that all of the 
EBL from sources with 850 /rm fluxes greater than 4 mJy 
is in the z = 1 — 3 range, where the identified sources 
of Table m lie. The fact that deep radio imaging only 
detects ~ 60% of the bright SCUBA sources suggests 
that the remaining 40% may be at redshifts greater than 
z ^ 3; thus, we may be overestimating this contribution. 
This result is plotted in Fig. d as the dashed line. The 
SFRD measured in this way is not independent of the 
SFRD measured in the core NIR sample, since the bright 
850 ^m sources may already be contained in the core NIR 
sample. Indeed, the directly determined value lies below 
the core NIR SFRD at these redshifts. This is because, 
while we are making the extreme assumption that all of 
the bright sources lie in this redshift range, only a small 
fraction of the 850 /rm EBL lies at these bright fluxes 
(about 7 Jy deg“^). 

Finally, the maximum completeness correction for the 
submillimeter SFRD in this redshift interval can be de¬ 
termined from the residual 850 fim EBL that is not iden¬ 
tified by the core NIR sample, assuming that it all lies 
at z ~ I — 3. Most likely only a part of this missing light 
will come from this redshift range and at least some may 
come from higher redshift sources where the counterparts 
at other wavelengths will be much fainter. This remain¬ 
ing light is approximately 7—21 Jy deg“^ in 850 /rm EBL, 
depending on the total EBL value used. The SFRD from 
the residual submillimeter EBL is shown by the rectan¬ 
gular area in Figure d and can be added to the SFRD 
from the core NIR sample to obtain the maximum pos¬ 


sible submillimeter SFRD at each redshift. 

The SFRD derived from the core NIR sample using 
the radio data is 0.09 Mq yr“^ Mpc“^ at z = I and is 
almost identical to that at z = 2 — 3. Even the maxi¬ 
mum completeness correction of 0.05 Mq yr“^ Mpc“^ to 
0.15 Mq yr“^ Mpc“^ can only result in a star-formation 
history that is rising slowly above z ~ 0.8. When in¬ 
tegrated through time, it appears that the majority of 
the star formation occurs around z = I. Thus, previous 
claims that the submillimeter SFRD is strongly peaked 
at z = 2 — 3 appear to be incorrect. 

5.3. Comparison with the UV-Determined Star 
Formation Rate Density 

In Figure d we compare the SFRD determined here 
with the SFRD determined from UV observations. Both 
are computed consistently with a Salpeter IMF extend¬ 
ing to O.I Mq. However, rather than plotting p versus 
redshift, as is usually done, we plot the quantity p x t 
versus the cosmic time, t. The advantage of this dis¬ 
play is that we can see more directly how many stars are 
formed at a given time. We use diamonds and triangles 
to show the direct star formation determinations from 
the UV light without any correction for extinction. We 
use solid squares to show our directly measured FIR star 
formation, and we use open squares to show our maxi¬ 
mal corrections for the missing EBL, if it is formed in 
the z = I — 3 redshift range. 

The FIR star formation is about a factor of 3 to 5 
higher than the uncorrected UV star formation, depend¬ 
ing on the correction applied for the missing 850 pm 
light (dotted curves). The lower value corresponds to the 
directly measured light, and the upper value to the max¬ 
imally corrected light. This is consistent with the usual 
dust corrections applied to the UV star formation rates, 
but it should be noted that the samples are somewhat 
disjoint, in that the blue star-forming galaxies (our tem¬ 
plate classes 6 and 7) contribute substantially to the UV 
star formation but not to the FIR star formation. Thus, 
the extinction corrections must be higher in the galaxies 
that produce the bulk of the FIR light and lower in the 
blue star-forming galaxies. 

We have also computed the total amount of stars 
formed as a function of time by combining the FIR and 
UV star-formation histories and integrating with respect 
to cosmic time. Above z = 4, where we do not have 
any information about the FIR light, we have assumed 
that the total star formation is five times the UV star 
formation, but the results at later times are quite insen¬ 
sitive to this assumption, since only a very small part 
of the total star formation occurs at these early times. 
The cumulative total is shown by the solid curves, where 
the lower curve corresponds to the directly measured FIR 
star formation and the upper curve to the maximally cor¬ 
rected FIR star formation. At z = 0, the curves match 
suprisingly well t o the local deter mination of the stel¬ 
lar mass density llGole et al.ll200lfl . which is shown by 
the large solid circle in Figure El The agreement is 
best for the lower estimates, but even for the maximally 
corrected case, the curve is only slightly high. For the 
directly measured case, half of the star formation occurs 
at z < 1.3, while in the maximally corrected case, this 
rises to z < 1.45. 
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Fig. 15.— Comparison of the star formation history vs. cosmic time from our FIR determinations with that from UV samples. We 
plot the star formation rate density computed for a Salpeter IMF extending to 0.1 Mq multiplied by the cosmic time, since this shows 
more clearly what fraction of the baryonic mass in stars is created at any time. The triangles and diamonds show the star formation 
that is directly seen at rest-frame UV wavelengths. At late times, we show the local Galaxy Evolution Explorer (GALEX) determination 
of Wyder et al. (2005) (large solid triangle), the GALEX determinations of Schiminovich et al. (2005) (small solid triangles), and the 
ground-based determinations of Wilson et al. (2002) (small solid diamonds), which are the most accurate measurements near z = 1. At 
intermediate times, we show the determinations of Steidel et al. (1999) (open diamonds), Bouwens, Broadhurst, Sz Illingworth (2003) (small 
open triangles), and Iwata et al. (2003) (solid upside-down triangle). The dashed curves show a parameterized fit. The dotted green curves 
show the total star formation that would be inferred if we were to apply reasonable extinction corrections to the UV light (upper curve = 
5, lower curve =3). The solid red squares show this paper’s FIR determinations, and the open red squares show our maximal corrections 
for incompleteness in the .2 = 1 — 3 range. There is broad agreement of the FIR determinations with the dust-corrected UV determinations, 
though the populations giving rise to the FIR light are somewhat disjoint from those giving rise to the UV light. The horizontal line shows 
the cosmic baryon density. The solid curves show the cumulative star formation. The solid black circle shows the present-day stellar baryon 
density estimated by Gole et al. (2001). 


6 . SUMMARY 

We have obtained accurate redshifts for the sources in 
the GOODS-N area using existing spectroscopic redshifts 
and improved photometric redshifts from NIR and MIR 
data. The radio-identified bright (> 2 mJy) SCUBA 
sources in this area are in the redshift range z ~ 1 — 3 
and have a median redshift of z = 2.5, consistent with 
previous radio and spectroscopic surveys. 

However, we used a stacking analysis to show that 
much of the 850 /im EBL is in fact traced by a NIR 
sample constructed from sources with fluxes greater than 
1.8 /iJy in either the H or 3.6 ^m bands. We showed that 
much of this light arises from galaxies with intermediate 
spectral types at z < 1.5. Thus, many of the fainter sub¬ 
millimeter sources that give rise to most of the 850 /rm 
EBL are at lower redshifts and lower luminosities than 
the bright submillimeter sources that are detected di¬ 
rectly. 

Finally, we used a stacking analysis to estimate the 
average 20 cm EBL produced by the unidentified or in¬ 
termediate spectral type galaxies in our NIR sample as 
a function of redshift, from which we determined the 
SERB. We found that this SERB evolves rapidly between 


z = 0 and z = 0.8, after which it becomes approximately 
flat. Using the submillimeter data directly, we then cal¬ 
culated a submillimeter based SERB at z ~ 1 — 3 which 
agrees closely with the radio based SERB. In addition, by 
assuming that all of the submillimeter EBL that is not 
accounted for by our NIR sample is also at these red¬ 
shifts, we put an upper bound on the SERB at z ~ 1 — 3. 
Even with this maximum completeness correction, we 
found consistency with a nearly flat or slowly rising ex¬ 
trapolation of the SERB from z ^ 1. We conclude that 
the majority of the star formation traced by the submil¬ 
limeter light comes from redshifts near one rather than 
at the higher redshifts that have been favored until now. 
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the Wisconsin Alumni Research Foundation, the Alfred 
P. Sloan Foundation, and the Bavid and Lucile Packard 
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